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We study bifurcations in a spatially extended nonlinear system representing population dynamics 
with the help of analytic calculations based on the time-independent Schrodinger equation for a 
quantum particle subjected to a uniform gravitational field. Despite the linear character of the 
Schrodinger equation, the result we obtain helps in the understanding of the onset of abrupt 
transitions leading to extinction of biological populations. The result is expressed in terms of Airy 
functions and sheds light on the behavior of bacteria in a Petri dish as well as of large animals such 
as rodents moving over a landscape. 
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I. INTRODUCTION 



The purpose of this paper is two-fold: to show how 
standard Schrodinger equation analysis familiar in a sim- 
ple physics situation may be used to understand an im- 
portant problem in biology, and to study, via such an 
analysis, the effect of spatial variation of resources on the 
ecologically relevant problem of extinction of populations 
of biological species. Extinction in biology is intimately 
related, as we shall show below, to abrupt transitions in 
physical systems, and is important whether the focus is 
on the dynamics of bacteria in a Petri dish [H, 0, H, 0, H, Q 
or of the movement of larger animals such as rodents in- 
volved in the spread of epidemics @, [1, E3] • 

Consider the logistic equation 



du(t) 
dt 



au{t) - bu 2 (t) 



(1) 



which may describe the time dependence of the density 
u of a spatially uniform population subject to a growth 
rate a and a resource competition term controlled by the 
parameter b. The steady state density undergoes a tran- 
scritical bifurcation at a — 0, as a is varied from negative 
to positive values. The physical meaning of the bifurca- 
tion is straightforward. Whereas for positive a the stable 
steady state value of u is a/b, for negative a the density 
vanishes in the steady state because the growth rate is 
then really a death rate. Introduction of space as a mere 
parameter X brings in no new behavior. 

A new aspect does enter, however, if spatial derivatives 
are introduced into Eq. ([T]) through, for instance, diffu- 
sion terms which represent the motion of the members of 
the species from one spatial location to another. It is this 
situation, represented by a Fisher equation [ll|, [l2|with 
spatially varying growth rate a(X) and diffusion constant 
D, 



du(X,t) 
dt 



a(X)u(X,t) - bu 2 {X,t)+D 



d 2 u{X,t) 
dX 2 



(2) 



that is the subject of the present paper. We study the 
1-d case for simplicity. As expected, the steady state 
u(X) has a tendency to follow the spatial variation of 
the growth rate a(X) but is modified by the effects of 
diffusion. For instance, if a(X) has the value a in a seg- 
ment of space, and the value — oo outside it, as shown 
first by Skellam (l3|, u(X) is given by the square of the 
elliptic cd function [5|, LL3[ . Thus u(X), which is maxi- 
mum at the center of the favorable region and drops off 
gradually towards the edges, is not proportional to a(X). 
The latter has the shape of a step function. 

Equation ([2]), which has been quite successful in the 
description of ecological dynamics [HI, [13], predicts in- 
teresting behavior in addition to the mere difference in 
shapes of a(X) and u{X). The maximum value u m of the 
density within the favorable region, i.e., the spatial seg- 
ment in which a(X) is positive, naturally depends on the 
length L of the segment, and decreases if L is decreased. 
What is striking is that there is a minimum value of L 
below which u m , and therefore the density throughout, 
vanishes identically. This minimum (critical) value can 
be obtained, except for a proportional constant equal to 
7T 2 , by equating the growth time l/o to the diffusion time 
L 2 /D: 



(3) 



The underlying physical idea is that the critical length of 
the favorable segment corresponds to the situation that 
the random walker (bacterium or animal) traverses the 
length of the segment diffusively in the time necessary 
for growth, and falls prey to the harsh conditions outside 
the segment. 

Our interest in the present paper is to ask how this 
well-known result [f| [T2, EH is modified by the simulta- 
neous introduction of two items of realism into the growth 
a(X): a gradual, rather than step-like, spatial variation 
of the growth rate at the border of the segment, and 
a non-infinite value for destruction outside the segment. 
The latter, by itself, has been analyzed earlier [lj] but 
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the former has not, to the best of our knowledge. We 
provide an analysis of the combined situation below on 
the basis of wave function calculations from the time- 
independent Schrodingcr equation describing a quantum 
particle in a uniform gravitational or electrostatic field. 

The analysis we present is applicable to at least two 
specific biological systems that have been studied re- 
cently: a bacterial population in a Petri dish, and rodents 
responsible for Hantavirus infection spread. In several 
experiments [l|, Q on the former system, the bacterial 
population is bathed in the presence of deadly ultravio- 
let radiation, except in a localized region covered with a 
mask that protects the bacteria. The growth rate a(X) 
is positive in the protected region and negative outside. 
In observations of the second system, there are patches 
of resources which make possible population growth and 
sustainance of rodents such as mice in parts of the land- 
scape but not beyond those patches. Earlier analyses 
have treated the patches or the masks to correspond to 
sharply changing resources. The real experimental sit- 
uation is naturally different. The ultraviolet light that 
falls on the mask seeps under the edges of the mask. As 
a consequence, the change in the growth rate is gradual 
rather than sudden. Similarly, the available water and 
food in the rodent patches change gradually in the land- 
scape. Transitions have been reported in experiments on 
the bacterial system Q following quantitative predictions 
made earlier |5(. Transitions also seem to be present in 
rodent experiments although the observations are more 
difficult to unravel. The observational relevance of our 
present investigation should be, thus, clear in both cases. 

The paper is set out as follows. A simple model that 
incorporates a gradual variation of resources in space is 
presented in section 2 and analytical formulae are derived 
in terms of Airy functions. Limits are studied in section 
3 and concluding remarks appear in section 4. 



II. MODEL AND THE RESULT 



The model we analyze is as given in Eq. ^ with the 
growth rate a(X) as in Fig.l. It equals the constant a for 
a spatial segment of length L around the origin X = 0, 
and drops gradually (linearly) to the constant negative 
value —a\ outside the segment, the spatial variation be- 
ing linear (constant slope). During the linear drop, the 
rate remains positive on either side of the segment for an 
additional distance R as shown. Using the suffixes 1, 2 
and 3 for the respective regions I, II and III in Fig.l, we 



have 

du r (X,t) 
dt 

du 2 (X,t) 
dt 

du 3 (X,t) 
dt 



= D 



D 



= D 



d 2 Ul {X,t) 

dX 2 
d 2 u 2 {X,t) 

dX 2 
bu 2 (X,t), 
d 2 u 3 (X,t) 
dX 2 



+ aux{X,t) - bu 2 (X,t), 



+ a 1 



2X - L 
2R 



u 2 {X,t) 



a lU3 (X,t)-bu 2 3 (X,t).(4) 



We do not attempt to solve either this nonlinear par- 
tial differential equation or the corresponding nonlinear 
ordinary differential equation for the steady state u(X) 
obtained by setting the time derivative equal to zero. 
Rather, we argue that, because our interest lies only in 
the extinction phenomenon, i.e., the vanishing of u(X), 
we can eliminate the bilinear terms of the density, and 
match solutions of the ensuing linear equation across the 
interfaces of the three regions I, II and III (see Fig. 1). 
This argument appears to have been given first by Lud- 
wig et al. [1J] in his investigations of the spatial pat- 
terning of the budworm. It has been used recently by 
the present authors to propose double-mask experiments 
for the bacterial problem [6|. The crux of the argument 
is that, at the transition, the densities vanish and there- 
fore terms of order higher than the first may be safely 
neglected even while the linear terms are compared to 
extract the extinction condition. The critical value of L 
signifying extinction which result for any positive b can 
therefore be calculated precisely by considering the lin- 
ear equation obtained by putting b equal to zero. The 
problem we are presented with here is simply of finding 
of the wavefunction for a time-independent Schrodinger 
equation of a quantum particle in a potential which is 
linear in X as in the case of a uniform gravitational or 
electrostatic field. 

Defining the dimensionless position x and the quan- 
tities I and r, all obtained by dividing A, L/2, and R 
respectively by the diffusion length yjD/a, we focus on 
the linear steady state counterpart of Eq. Q, with the 
notation £ = y/aija: 

\ + Ul (x ) =0, 

ax 



d 2 u 2 (x) 



dx 2 

, ?u 3 (x) = 



i - ) u 2 (x) = 0, 



d 2 u 3 (x) o 



dx 2 



(•5) 



The growth rate plays the role of the potential in a 1- 
dimensional quantum mechanical problem in a constant 
force field. The first of these equations is solved trivially 
in terms of trigonometric functions, the third in terms of 
exponentials, and the second in terms of Airy functions. 

Specifically, using the fact that the system has sym- 
metry about x = 0, the solution in the three regions is 



3 



a(X)/ 





where Ai and Bi are Airy functions and the A's and 
B 2 are constants to be determined from the boundary 
conditions. 



FIG. 1: The spatial dependence of the growth rate a(X) show- 
ing its value as the constant a in a segment of length L joined 
through a linear decrease to the value —01 on both sides. The 
width of the transition region is R on either side. 



given by 



U\(x) 
u 2 (x) 



Aicos{x), 



= A 2 Ai 



B 2 Bi 



u 3 (x) = A 3 e ix , 



r 2/3 h 

_ r 2/3 , , 



x-l 
r 

x — I 



(6) 



Matching the logarithmic derivative of the solutions 
across the interfaces and introducing rj = r 1 / 3 for nota- 
tional convenience, we arrive at the following condition 
for extinction. 



tan / = 




Ai'{^rf)Bi'{-if) - Ai'(-r) 2 )Bi'(t 2 v 2 )+Zr) [Ai(t 2 r) 2 )Bi'(-r?) - Bi(frj 2 )Ai' \-7 1 2 )} 
Aii-^Bi'i^r] 2 ) - Ai'{ev 2 )Bi{-V 2 ) + £77 [Ai(-?7 2 )Bi(£V) - Ai^ 2 V 2 )Bi{-r] 2 )} 



r 



(7) 



Primes denote derivatives 
length I 



The normalized extinction 



L/y/D/a is the arctangent of the right hand 
side if the latter is positive, and is zero otherwise. 

Equation ([7]) is the central result of this paper. The 
crucial parameters in the expression are 77 which mea- 
sures the transition region relative to the diffusion length 
and £ which measures the death rate outside, relative to 
the growth rate inside, the segment. As is clear from 
the expression, the second parameter can also be chosen 
naturally as the product £?y . 



III. LIMITING CASES 

If the transition region vanishes (r — 0), we take the 
limit 77 — > in Eq. ([7]) and recover the Ludwig formula 

GJ 



L = Ld = 2\ — arctan w — . 

V a \ a 



(8) 



In the general case r ^ 0, the extinction value normal- 
ized to this Ludwig value is given by plotting the analytic 
expression in Eq. ©■ Fig. 2 shows, in the form of the 



solid line in the main figure, this dependence of L/Ld on 
r, the normalized extent of the transition region, for a 
specific value of £, viz., 3. The solid circles are obtained 
via the numerical solution of the full nonlinear equation 
(j4]), for a nonzero b — 1 in appropriate units. Note that 
Eq. (HJ) cannot be solved analytically, and that the com- 
plete agreement between the numerical results from the 
nonlinear equation (whatever the value of b) and the an- 
alytic result from the linear extinction analysis (that has 
no b) vindicates our theoretical procedure. 

Various limiting cases may be studied from our Eq. 
through ascending and asymptotic expansions of the 
Airy functions. An example of the former is (T5 , [TtJ 



Ai(x) = ci/(x) - c 2 g(x), 
Bi{x) 



V3 



where, 



cif(x) + c 2 g(x), 



1 3 1A 6 
-.x 6 -I -x b 



(9) 
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by the ascending series approximation 



denoted by a dashed line, with c = — Ai(0) / Ai' (0) = 
1.372. Important to note is that the exact L becomes 
zero at a finite value of the transition region extent (at 
r = 1.028). This follows from the existence of a zero of 
Ai'{— r 2 / 3 ) at finite r; the physical significance is com- 
mented on below. 



FIG. 2: Excellent agreement between our Airy extinction pre- 
diction formula (0 based on the linear theory (solid line) and 
the numerical solution of the full nonlinear equation (circles). 
Plotted in the main figure is the extinction length normalized 
to the Ludwig value, L/Ld, versus the extent of the transition 
region r for £ = 3.0. The inset shows steady state density pro- 
files corresponding to three different values of L larger than 
the critical, for £ = 3.0 an d r = 0.4. Distance x is expressed in 
units of diffusion length y/D/a while the steady state density 
u(x) is expressed in units of a/6. The density decreases as the 
patch length decreases and vanishes at the critical value, at 
which L/L d = 0.398. 




FIG. 3: Generalization of the Skellam result © in the pres- 
ence of a transition region of extent R. We plot L/L 3 as a 
function of r = R/^JTTJa. The solid line is the exact result 
(Eq. UOp and the dashed line the approximation retaining the 
lowest power of r 1/3 . Note that L = for r > 1.028. 



and ci = Ai(0), c 2 = Ai'(0). 

One of the goals of the present paper is the gener- 
alization of the Skellam result ((3]) in the presence of a 
transition region. For this, we let £ — > oo in (JT]) and 
obtain 



L = 2\l — arctan 
a 



1 Ai'( 



_ r 2/3 



r i/3 Ai(-r 2 / 3 ) 



(10) 
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FIG. 4: Our analytic predictions for the dependence of the 
extinction length L (in units of ^/D/a) on the transition re- 
gion extent in (a) and against the ration of the death to the 
growth rate in (b), plotted from the Airy function formula 
Eq. (0. Note the threshold values of 1/r and £ below which 
extinction does not occur. 



The extinction length L given by the exact expression 
([7|) is plotted in Fig. 4 against (a) the reciprocal of the 
transition region extent, 1/r, for different values of £ as 
shown, and (b) the square root of the ratio of the death 
rate to the growth rate, £ for different values of r. The 
extinction length becomes zero for large enough r in (a) 
and small enough £ in (b). This means that there is no 
extinction for r > r c and £ < £ c where r c and £ c are the 
threshold values. These threshold values are plotted in 
Fig. 5 directly from the graphs. We have verified that the 
curves coincide precisely with those obtained by putting 
the numerator of Eq. ([7]) equal to zero. 




FIG. 5: Plots of threshold values £ c and r c obtained graphi- 
cally. They are found to agree with the prediction of Eq. ([7} 
obtained by putting its numerator equal to zero. 




We divide Eq. §TU§ by the Skellam result ©, and plot 
L/L s MS i\ function of the transition region extent (i.e., 
r = R/ y/D/a) in Fig. 3. The exact expression is the 
solid line. It is approximated reasonably well for small r 



This vanishing of the extinction length L, signify- 
ing that the extinction phenomenon disappears for large 
enough transition region or small enough death rate, 
stems from the definition of L as the top segment of the 
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a(x) trapezoid in Fig. 1. Thus, Fig. 6 shows with this 
segment of zero extent (L = 0) that the favorable patch 
is still wide enough to make possible for the walkers to 
survive in the protected patch because of the transition 
region. 



a(x) 
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a 


"I / 


o \ 





X 



FIG. 6: Situation where L = 0, there being no extinction 
because the transition region itself provides a large enough 
favorable patch for survival. 



IV. CONCLUSION 

Our focus in this paper has been to study the combined 
effect of a gradual change in resources at the borders of a 
patch, and of finite death rates outside the patch, on the 
extinction transition of bacteria in a Petri dish, rodents 
in a landscape, or of similar random walkers that can 
move over an inhomogeneous spatial region. Although 



the system is essentially nonlinear, and the extinction 
transition is a consequence of the nonlinearity (coupled 
to the inhomogeneity), we have been successful in ob- 
taining exact solutions through an analytic procedure. 
The procedure draws on an analogy to the solution of a 
time-independent Schrodinger equation and results in an 
expression for the critical length of the patch in terms of 
Airy functions and their derivatives. 

Our central result is Eq. ([7]). One of the crucial pa- 
rameters is r\ which is r 1 / 3 , where r is the length of the 
transition region in units of the diffusion length y/D/a. 
The appearance of the latter is natural since the under- 
lying process is the movement of the random walker (the 
diffusing entity) from within the favorable patch to the 
unfavorable space outside within the growth time. The 
emergence of the third power is a consequence of the Airy 
function which arises from the linear dependence of the 
gradual decrease of the resources. 

The other quantity of importance is the ratio of the 
death rate outsid e the patch to the growth rate inside, 
specifically £ = yj a±/a. This second parameter can also 

be taken to be the combination £77 = -\/r 2 / 3 ai/a. We 
have recovered the Skellam result Eq. ([3]) and given its 
generalization in Eq. (JTDJ). We have also recovered the 
Ludwig result and provided its generalization as well. 
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